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Abstract 

We compare the explicitly correlated Hylleraas and exponential basis sets in the evaluations of ground 
state of Li and Be"''. Calculations with Hylleraas functions are numerically stable and can be performed 
with the large number of basis functions. Our results for ground state energies —7.478 060 323 910 10(32), 
— 14.324 763176 79043(22) of Li and Be+ correspondingly, are the most accurate to date. When small 
basis set is considered, explicitly correlated exponential functions are much more effective. With only 128 
functions we obtained about 10^^ relative accuracy, but the severe numerical instabilities make this basis 
costly in the evaluation. 
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I. INTRODUCTION 



In order to accurately calculate energy levels of light atomic systems, not only nonrelativistic 
energies, but also relativistic and QED corrections have to be obtained with the high precision. In 
the NRQED approach all corrections are obtained perturbatively, in powers of the fine structure 
constant a. Each term of this expansion is expressed as the expectation value of some effective 
Hamiltonian with the nonrelativistic wave function. Similarly, corrections due to the finite nuclear 
mass and its size can all be included perturbatively. This however requires the accurate represen- 
tation of the nonrelativistic wave function. 

The wave function of the ground and excited states can be obtained on the base of the Ritz 
variational principle. The accuracy of the upper bound for energy mainly depends on the basis set 
of trial functions and effectiveness of the optimization routine. There are not so many possible 
choices of basis functions, knowing that electron correlations have to be accurately accounted for. 
The most serious problem in development of explicitly correlated methods is difficulty in accurate 
calculations of integrals appearing in Hamiltonian matrix elements, and the complexity of these 
integrals grows with the increasing number of correlated electrons. 

The most often in use are correlated Gaussian functions which have been applied so far to sys- 
tems including up to six-electrons, and the most accurate results in comparison to other methods, 
have be obtained for Be atom . Relatively simple integrals and possible generalization to 

systems with higher number of electrons is the main advantage of Gaussian functions. However, 
these functions have improper short-distance (Kato cusps) and long-rage behavior. As a result, 
the convergence of the variational procedure is not very fast. Quality of the globally optimized 
trial functions, even in a few thousand basis set is often insufficient for calculations of relativistic 
effects beyond the leading order. In particular, we observe poor convergence of matrix elements 
with singular operators i.e. Dirac b. 

Until now, the most accurate nonrelativistic wave function for lithium-like atomic systems were 
computed in Hylleraas basis by King in [4], by Yan and Drake in [^Sfl and by present authors in [^. 
The Hylleraas function for the three-electron system is of the form 

with nonnegative integer values of Ui. Although, algorithms for integrals with these functions are 
computationally demanding, the correct long and short-range asymptotic and possibility to use a 
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large basis set of functions (~ 10000) with small number of variational parameter (~15) allows 
one to achieve high accuracy. In a recent series of papers we formulated the analytical method for 
calculations of Hylleraas integrals with the help of recursion relations In this work we tuned up 
the optimization routine compared to our former work [6]. As a result, we significantly improved 
nonrelativistic energies as compared to the previously published ones in [SilSQ and achieved about 
10^^'' precision. 

Even better precision can in principle be achieved with the explicitly correlated exponential 
function. In 1987 Fromm and Hill obtained the closed analytical formula for the related four- 
particle integral 

/(Pri f (Pr2 f (Pr^ Q-wiri-W2r2-W3r3-uir23-U2ri3-U3ri2 
~A / 1 / 1 ' 

An J An J An ^23 rgi ri2 ri ra rg 

reducing the problem to the evaluation of multivalued dilogarithmic functions of complex argu- 
ments [9] . Their formula could be differentiated with respect to the Wa and Ua to introduce pre- 
exponential powers of the and Tab, thus to generate the class of integrals needed for evaluation 
of Hamiltonian and overlap matrix elements. The Fromm-Hill formula was modified later by 

o 

Harris eliminating the necessity of branch tracking on the complex plane [HOD . Zotev and Re- 



bane presented their method for integrals with an extension to complex exponentials llllh . They 



demonstrated fast convergence even in small bases and high potential of this method in variational 



calculations of four-body systems [|12|] . Recently, Guevara et al. [13] have been able to optimize 
the correlated exponential function including linear terms in inter-particle distances by the six- 
dimensional numerical integration and obtained nonrelativistic energy with the relative precision 
of about 10^1 

Effectiveness of correlated exponential functions gives opportunity to reduce significantly the 
size of the basis set as compared to Gaussian and Hylleraas functions. However, the evaluation 
of corresponding integrals is the most time consuming part of the variational method. This fact 
suggests to use rather short basis with carefully optimized parameters. In this work these integrals 
are calculated as folows. The master integral qq in Eq. ^ is calculated using Harris formula 
[fioll . Integrals with higher powers of inter-particle distances, are obtained using recursion rela- 



tions, which are derived from the differential equation (|T8l) . As a demonstration of this method, we 
performed numerical calculations of the nonrelativistic energy and of Dirac-(5 for the ground state 
of Li and Be"*". With 128 well optimized correlated exponential functions with real parameters we 
have obtained nonrelativistic energies with relative precision of about 10^^. This precision is not 



impressive in comparison to the value extrapolated from 13944 Hylleraas functions. However, the 
result for lithium is comparable to six times bigger set of Hylleraas functions or 1500 optimized 
Gaussians. The highly accurate wave function in a small basis set gives a flexibility in devel- 
opment of numerical methods for evaluation of more complicated integrals. It is expected to be 
especially valuable for evaluation of matrix elements of m operators, which involves integrals 
very difficult to deal with Hylleraas functions. 



II. NONRELATIVISTIC WAVE FUNCTION 

The ground state wave function ^ is represented as a linear combination of ip, the antisym- 
metrized product of the spatial functions and the spin function % 

ijj = A[(j)(fi,f2,f3)x] , (3) 

X = a(l)/?(2)a(3)-/?(l)a(2)a(3). (4) 
In the case of correlated exponential functions, 0(ri, r2, rs) is 

and we assume that ctj, Pi are real numbers. These nonlinear parameters are subject of additional 
conditions. Namely, when one of the electrons goes to infinity, the wave function shall decay 



exponentially sufficiently fast, so for example ai + + /^s > V^^E~, where E'ion is the ionization 
energy. 

The expansion coefficients and nonlinear parameter are obtained by minimization of energy 
with the Hamiltonian H 

H ^ T + V, (6) 

1 „ 1 ^ J, 1 ttO 

a=l a=l a>b=l 

where Z e is the nuclear charge and atomic units are used elsewhere. After elimination of spin 
variables, the matrix element of H can be expressed as 

= (2 0^(1, 2, 3) + 2 <^^(2, 1, 3) - 0^(3, 1, 2) - 0^(2, 3, 1) - 0^(1, 3, 2) 

-(/.^(3,2,l)|i/|(^^(l,2,3)). (8) 
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The individual matrix element {(f)^\H\(j)^) is represented as a linear combination of 34 Slater 
integrals defined as 

J An J An J An 

^"1-1 ^"3-1 »,n4-l ^"5-1 »,n6-l fQ\ 

'23 '31 '12 '1 '2 '3 ) 

where are nonnegative integers and Wa = + , = /Jf' + P^. The number of necessary 



integrals for the matrix element of H can be significantly reduced. Rebane and Zotev [|l4h derived 
the formula which includes only seven integrals: the overlap integral and six Coulomb 

integrals (0^ |r~^ 1 0^) , which we have found very useful. It reduces significantly the computational 
costs in most of cases except for small Wa, Ua, where it becomes numerically unstable. In this case 
we use the numerically stable standard form of the kinetic energy operator obtained by direct 
differentiation of the left and the right wave function over the electron coordinates. 

III. CALCULATION OF SLATER INTEGRALS 
A. Integration by parts method 

The evaluation method of g{ni, n2, n^, n^, n^, uq) in Eq. ^ is based on the integration by parts 
identities, which are widely used for the analytical calculation of Feynman diagrams [15]. Let us 
consider the following integral in the momentum space 

G{mi,m2,m3; 1714,1715, me) = ^ j dPki j d^ki J d^ ks {kj + uj)~"^^ {kj + ul)'""^ 

{kl + ul)-"^^ {kl, + wlr-^ {kl + wl)-^'" {kl + wlr-tlO) 

which is related to g function by go = g{0, 0, 0, 0, 0, 0) = G{1, 1, 1, 1, 1, 1). There are 9 corre- 
sponding integration by parts identities 



d 



,2\— mi 



= id(z,i) = J d'ki J d^k2 j d-'ks {kf + uf 

{kl + Ulr^^ (kl + Ulr^^ikl, + (fc?3 + u;2)-m. (^2^ ^ ^2)-me 



(11) 



where i,j = 1,2,3. The reduction of the scalar products from the numerator leads to the relations 
between functions G of different arguments. These identities group naturally into three sets with 
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respect to j. For example for j = 3 and mj = 1 we have the following system of three equations 



G(0, 1, 1, 1,2, 1) - ^(0, 1,2, 1,1,1) + G'(1,0, 1,2, 1,1)-G'(1, 1,0, 1,2,1) 
-G(l, 1,0,2, 1, 1) + G(l, 1, 1, 2, 0,1)-G'(1, 1,1,2, 1,0) + ^(1, 1,2, 1,0,1) 



+G(1, 1, 1, 1, 2, 1) i-uj + 4- wl) + G(l, 1, 2, 1, 1, 1) (uj + u 



Wo 



+G(1, 1, 1, 2, 1, 1) i-ul + 4- wl + wl). 

G(0, 1, 1, 1,2, 1) + G(1,0, 1,2, 1,1)-G'(1,0,2,1,1,1)-G'(1, 1,0, 1,2,1) 
-G(l, 1,0,2, 1, 1) + ^(1, 1, 1,0,2, 1)-G'(1, 1,1, 1,2, 0) + G'(l, 1,2, 0,1,1) 
+G(1, 1,1,2, 1, 1) {-ul + ul- wl) + ^(1, 1, 2, 1, 1, 1) {ul + ul- wl) 
+G{1, 1,1,1, 2, 1) {-u\ + ul-wl + wl), 

^(0, 1, 1, 1,2, 1) + G(1,0, 1,2, 1,1)-G(1, 1,0, 1,2, 1)-G'(1, 1,0,2, 1,1) 
-G(l, 1, 1, 1, 1, 1) + 2 ^(1, 1, 2, 1, 1, 1) ul + G(l, 1, 1, 2, 1, 1) {-ul + ul + wl) 
+G{lA,l,l,2,l){-ul + ul + wl), 



(12) 



Whenever mj = 0, G becomes a known two-electron integral F as defined in Appendix A. For 
example 



G(0, 1, 1; 1, 1, 1) = T{-l,Q,-l;w2 + W'i,Wi,U2 + us) 

- \u(l - ^2 + ^3 + W2 + W3 

2Wi\_ \ U2 + U'i + Wl 

^1 f Uh + W2 + W^\ ^ 

2 V -Uo + M3 + -Wl / 6 



)+u(i 



U2 + U2, + W2 + Ws\ 
Wi+W2 + W'i J 



U2 + U2, + Wl 

We solve the system of equation (fT2l) . for example against ^(1, 1, 1; 2, 1, 1), and obtain 

G{1, 1, 1; 1, 1, 1) - 2 a ^(1, 1, 1; 2, 1, 1) + P = , 



(13) 



2 dwi 

where cr is a polynomial 



(14) 



a 



I I I X III, 111, III, 1 I I I , I 

U^ U2 W^ + U2 U3 W^ + Ul ttg W2 + W2 W^ + Ul W-i [Ul + Wl 



Uo — Un — Wn — Wn 



,1111,1 1 1 1 1\ , 1 1 I 1 , 1 1 1 1 1\ fi i:\ 

+U2 W2 {U2 +W2-Ui-U.s-Wi- W^) + % W^ (% + W.s-U2-Ui-Wi- Wg) , (15) 
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and P is a the sum of two-electron integrals T 



P 



-Ui Wi [{Ui + W2)'^ - ul] r(0, 0, -1; Ml + W2, Us, U2 + Wi) 

-ui wi [{ui + Ms)^ - wl] r(0, 0, -1; Ul + M3, W2, Wi + W3) 

+ [Ui w-i + U2W2 — % W3 + Wi W2 [Ul + U2 — w^)\ r(0, 0, —1; Wi + W2, W3, Ul + M2) 

+ [ul wf - ulwl + ulwl + Wi {ul + ^3 - wl)] r(0, 0, -1; Wi + W^, W2,Ui + U3) 

-[u2 {u2 + wi) {u\ + u\ - wl) - u\ {u\ + u\ - wl)] 1(0, 0, -1; U2 + Wi, Us, Ul + W2) 
-[u3 (-U3 + Wl) {ul + u\- wl) - u\ {u\ + u\ - wl)] r(0, 0, -1; M3 + Wi, U2, Ui + 103) 

+Wi [W2 {u\ - u\ + wl) + Wz {u\ +wl- ul)] r(0, 0, -1; W2 + W^, Wi,U2 + Us) 

+wi [u2 {uj - wl + ul) + Us {ul + ul - wl)] r(0, 0, -1;U2 + Wi, W2 + W3) . (16) 



Since 



G(l,l,l;2,l,l) 
Eq. (fT4)) takes the form of a differential equation 



1 dgo 



2wi dwi 



dgo I da , TD n 
owi 2 owi 



(17) 



(18) 



or 



d 

v^^(v^i?o) + ^ = . 



(19) 



Analogous differential equation with respect to other parameters Wi and Ui can be obtain by ap- 
propriate permutation of arguments, using the tetrahedral symmetry of the function qq. This dif- 
ferential equation has been previously derived in Ref. 



B. Calculation of 



go was obtained in analytical form by Fromm and Hill in 



0] 



in terms of combination of multi- 



valued dilogarithmic function of complex arguments. Their formula was later simplified by Harris 
lioll . who was able to eliminate the ambiguity of choosing the right branch of dilogarithmic func- 
tion. In this work we use directly his formulae and allowed ourselves to verify its correctness. For 
this we used the solution of the differential equation in terms of one-dimensional integral. Namely, 
for a > we find 

g, = ^( r dw'i^^t^ + goy^l ), (20) 
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where 

9o 



sgn(wi) 



toi=oo 



TT 



— + - In^ 
6 2 



^41 + ^3 + W2 
Ui + U2+ Ws 



Li2 1 



U2 + U3 + W2 + W3 

U1+U3 + W2 



+ 



Li2(^] 



U2 + U3 + W2 + W3 



(21) 



U1 + U2 + W3 

The above integration over wi is performed numerically using adapted Gaussian points for the 
logarithmic singularity at wi = 00, see Appendix of [js]. 
For (T < we find 



9o 



-a 



Jill 



(22) 



where 



a 



w\=w\ 



0. 



(23) 



This integral is performed numerically using Gauss-Legendre quadrature in variable t 
y/wi — wi. In the simplest case when a = 0, go can be readily obtained from Eq. (fTSi) 



90 



-2P 



da 
dwi 



(24) 



In almost all the cases, we achieved 28 digits accuracy using quadruple precision arithmetic with 
about 100 integration points. 



C. Recurrence scheme 

Since the direct evaluation of g{ni, n2, n^, n^, 77,5, uq) in Eq. Q is very time consuming, it is 
desirable to derive recurrence relations permitting integrals of larger index values to be expressed 
in terms of those with smaller indices. From differential equation (fT9l) we can deduce much more 
than only integral representation for qq. We notice that 

9{ni,n2,n3, Hi, 115,716) = (-l)"'^'"'^"'^;^ • • • fi'o- (25) 

Analogously, we introduce (T(ni, 722, ris, ^4, 725, ng) and P(ni, n2, na, 714, ns, ng) derived form a 
and P respectively. If a 7^ then equation (fTSi) takes the form 

^ 0, 0, 0, 0, 0) ^(0, 0, 0, 0, 0, 0) + a(0, 0, 0, 0, 0, 0) 9(1, 0, 0, 0, 0, 0) = P(0, 0, 0, 0, 0, 0) . (26) 

Clearly this algebraic equation can be used to obtain 51(1,0,0,0,0,0) once 5^(0,0,0,0,0,0) is 
evaluated from the direct Ref. [Q] or integral (I20I22|) formulae. Now, we differentiate equation 
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261) ni — 1, 722, 713, ^4, ^5, timcs over wi,W2, ws, ui, U2, respectively 



Y] ( .M ••( (j{ni-ii,...,ne-i6)g{ii,--,i6) = P{ni-l,n2,n3,ni,nri,nG),(27) 

where we introduced a Newton-like notation 

?7,\ 1 /n\ ^ /ri\ /n — 1\ _^ ~ ^\ ^28) 

0/1/2 2' V^A/2 ' vvi/2 V ^ A/2 v-i/1/2' 

The above formula allows to express the integral g{ni, .., hq) with nonzero ni through (7-integrals 
with smaller index values. The expression for a{ni, ?t,2, ^3, 77,4, 71,5, tj-q) can be explicitely generated 
as derivatives of the polynomial a, since they become zero for large values of indices rii. P 
has a simple structure in terms of two-electron integrals T multiplied by a simple polynomial. 
Derivatives of these polynomials can be calculated explicitly. For T we use the recurrence scheme 
proposed by Korobov in IitI. 

Similar recurrence relations can be obtained from the differential equation like that in Eq. (fTSl ). 
but with respect to a different variable. We use them for the missing integrals with rz-i = in the 
above wi scheme, thus completing the algorithm for all (7-integrals starting from the master one 
go. We use them also to check the numerical stability of the recurrence scheme, as (^(1, 1,1, 1,1,1) 
can be obtained from the differential equation in any of these nonlinear parameters. As the result 
of this checking, we found out, that these recursions become unstable for small values of a in Eq. 
(fTSl ) and as a remedy we used higher precision arithmetics in this particular region. 

Recently, Harris obtained a family of recurrence formulas which enable construction of corre- 



lated exponential integrals with arbitrary pre-exponential powers of inter-particle distances [II8II . 
In comparison to them, our recurrences are not equivalent. Harris's recurrences in the denominator 
involve additional powers of Ui and thus may become numerically unstable in the limit of small 
Ui. This however, requires numerical verification. 

IV. OPTIMIZATION AND RESULTS 
A. Hylleraas basis set 

In Table I we present results obtained with Hylleraas functions for ground states of Li and 
Be+, as they are much more accurate than previous ones in [|5i|6|]. In comparison to these former 
works, we used slightly different division into 5 sectors with its own set of nonlinear parameters 



as proposed in Ref. [|5D, and enhanced the optimization process by replacement of the minimiza- 
tion routine with CG Polak-Ribberie I19I1 with modifications of the line search algorithm jlQ ]. 
In Ref. [0] we performed optimization in quadruple precision arithmetics. Here we observe that 
this precision is sufficient for determination of the nonrelativistic energy, but it is at the edge of 
numerical stability for analytical calculation of gradients in a basis set corresponding to = 
max(^. Hi) = 10. Therefore, in this work we used sextuple precision arithmetics for the whole 
calculation. Obviously, optimization process in higher precision arithmetics takes more time, in 
this case it is about 5 times longer, but the accuracy is improved by at least an order of magni- 
tude. The results presented in Table H] are better than the former ones in 50 percent bigger basis 
set. Especially important is the numerical result for maximum set of 13944 carefully optimized 
functions, as this guarantees good quality of extrapolation to 00 and estimation of an uncertainty. 



TABLE I: Ground state nonrelativistic energies for the ground state 
with Hylleraas functions with comparison to earlier results including 



of Li and Be+ for various basis length 
; correlated Gaussian functions. 



No. of terms 


£;(Li) 


^(Be+) 


2625 


-7.478 060 323 570 509 


-14.324763 176 517134 


4172 


-7.478 060 323 845 785 


-14.324763 176 746 865 


6412 


-7.478 060 323 898 268 


-14.324763 176 783 625 


9576 


-7.478 060 323 907 743 


-14.324763 176 789144 


13944 


-7.478 060 323 909 560 


-14.324763 176 790150 


00 


-7.478 060 323 91010(32) 


-14.324763 176 79043(22) 


9577^ 


-7.478 060 323 892 4 


-14.324763 176 766 8 


00^ 


-7.478 060323 906(8) 


-14.324763 176 784(11) 


10000^= 


-7.478 060 323 81 





8000'^ 
16764^ 



-7.478 060 3234519 



-14.324 763 1764 



a - Ref. 



b - Ref. m, c - Ref. tZE, d - Ref 



e - Ref. i2M 
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B. Correlated exponential basis set 

We optimized the correlated exponential basis set incrementally starting from 1 up to 128 func- 
tions as shown in Tables HI] and [nil At the starting point, the bigger basis was composed of 

TABLE II: Nonrelativistic energies and Dirac-(5 expectation values for the ground state of Li compared to 
results in Hylleraas basis 



N 


E(U) 


5E/E 






1 


-7.453 907 382 


3.210- 


-3 


13.631 327 


0.614 377 


2 


-7.465 318 352 


1.710- 


-3 


13.163 649 


0.617 596 


4 


-7.476 009761 


2.710" 




13.691 905 


0.586 670 


8 


-7.476 936 884 


1.510- 


-4 


13.773 519 


0.576 457 


16 


-7.478 052 680 


1.010- 


6 


13.840924 


0.545 361 


32 


-7.478 059401 


1.210- 


-7 


13.841641 


0.544 671 


64 


-7.478 060050 


3.710- 


-8 


13.842 162 


0.544 526 


96 


-7.478 060272 


7.010- 


9 


13.842 641 


0.544 391 


128 


-7.478 060301 


3.110" 


'9 


13.842 618 


0.544 368 


Hyll. oo 


-7.478 060 323 9 






13.842 610 8 


! 0.544324 6 



previously optimized smaller basis and functions with randomly chosen nonlinear parameters un- 
der constraints resulting from interparticle separation conditions. Due to the presence of many 
nonlinear parameters, each function has its own set of 6 parameters, the optimization process was 
divided into steps. In a single step nonlinear parameters of only one function were optimized using 
Powell method without gradient. In one cycle all functions were optimized separately. For small 
basis several cycles were needed to achieve convergence at the 9th digit after the decimal point, 
and for larger set of functions number of cycles increases. Implementation is done in Fortran 95 
in the quadruple precision arithmetics. In the region of typical values of Wa and we observe 
very good numerical stability of recurrence relations. However, in some particular cases during 
the minimization process, where a in Eq. (fT5l) becomes small and changes its sign, the sextuple 
precision arithmetics was needed, as the recurrence relations lose numerical precision. The region 
of small a is numerically unstable and we have not found yet an alternative way of evaluation of 
g{ni,n2, n^, n^, n^, uq) functions, by avoiding the presence of a in the denominator. This would 
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TABLE III: Nonrelativistic energies and Dirac-(5 expectation values for the ground state of Be+ compared 
to results in Hylleraas basis 



N 


T— 1 /T-) -|-\ 


SE/E S{ra) 


o{rab) 


1 


-14.269 015 274 


3.9 10-^ 34.584 174 


1.726084 


2 


-14.319 868 303 


3.410-'' 34.818 880 


1.722 376 


4 


-14.324 097 014 


4.710-5 35.163 138 


1.598 315 


8 


-14.324 646 319 


8.210-6 35.082068 


1.589484 


16 


-14.324 730041 


2.310-6 35.118 928 


1.583 949 


32 


-14.324760432 


1.910-'^ 35.109 851 


1.582 886 


64 


-14.324762726 


3.110-^ 35.102 872 


1.581 131 


96 


-14.324 763 106 


4.910-9 35.105 550 


1.580752 


128 


-14.324763 141 


2.510-9 35.105 342 


1.580583 



Hyll. oo -14.324763 176 8 35.105 055 7 1.580538 6 



be necessary for larger basis set and for states with the higher angular momentum. The quadruple 
precision arithmetics for the maximum basis of 128 functions guarantees high quality of the total 
wave function and the energy. We observe by comparison with Hylleraas results, that the relative 
accuracy of about 10^^ is achieved for energies, and about 5-6 significant digits for wave functions 
as indicated by the Dirac 6 expectation values. 



V. SUMMARY 



We have performed accurate calculations of the ground state energy and the wave function of 
Li and Be+ using explicitly correlated Hylleraas and exponential basis sets. Obtained results with 
Hylleraas basis are the most accurate to date, due to the use of large number of functions and 
efficient optimization. Results with correlated exponential functions are much less accurate, but 
they are the most efficient for the limited number of functions. The relative accuracy of about IQ-^ 
for the nonrelativistic energy of the ground state of Li and Be"^ with only 128 functions confirms 
high effectiveness of this basis. Compared to both the Hylleraas and the Gaussians functions, it 
allows to reduce significantly the size of basis set. Using the computational method based on 
recurrence relations, we are able for the first time to perform optimization process with as much 
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as 128 correlated exppnential functions and even more, if numerical instabilities for small a are 
eliminated, probably by a different type of recurrences. 

Our primary motivation for developing explicitly correlated exponential basis set is the efficient 
representation of the wave function in a small number of basis functions. We aim to apply them for 
numerical calculations of expectation values of operators corresponding to higher order relativistic 
and QED effects. They involve integrals with quadratic inverse powers of at least two interparticle 
distances. That kind of integrals are very complicated in the evaluation in the Hylleraas basis set 
and have not yet been worked out by the recursion method of the authors. However, there is a know 
algorithm by King [4], but his method is much too slow for a large scale computation. In the case of 
Slater integrals the problem would even much more complicated, but we think, one shall be able 
to perform this class of integrals numerically. Equiped with the large and accurately optimized 



Hylleraas basis [|24l1 . and with the short and flexible correlated exponential basis functions, we 
are aiming to determine m and m ol' effects in the hyperfine and fine structure of lithium-like 
systems. 
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APPENDIX A: TWO-ELECTRON INTEGRALS 

The two-electron integral T is defined by 



r(ni,n2,n3,Q!,/3,7) = 




ari-l3r2-'yri2 ^"i-l ^"-2-1 ^na-l 
'1 '2 '12 



(Al) 



An J 47r 



This integral takes very simple form when all = 



r(0,0,0,a,/3,7) 



1 



(A2) 



{a + P){a + 7) (/3 + 7) • 
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The explicit form for Ui > can be obtained by differentiation with respect to the corresponding 
nonlinear parameter, the result for negative rii is obtained by an integration, for example 

r(0, 0, -1, a, 7) = 7 ^ In f ^1 • (A3) 

[a -(3) [a + (3) \-i + 13 J 

For the actual evaluation of V we use compact reccurence relations from the work of Korobov 

0. 
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